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ABSTRACT 


The  general  acceptance/rejection  algorithm  for  generating  random 
values  on  a computer  is  specialized  for  distribution  tails,  using  an 
exponential  majorizing  function  and  a linear  minorizing  function.  Speci- 
fic algorithms  are  given  for  the  normal,  gamma,  Weibull  and  beta  distri- 
butions. While  the  algorithms  can  be  used  alone,  it  is  anticipated  that 
their  major  v_lue  will  be  to  serve  as  components  of  algorithms  for  complete 
distributions. 
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1.  INTRODUCTION 


Distribution  tails  are  treated  differently  than  the  rest  of  the  distribu- 
tion in  many  random  variate  generation  algorithms.  For  example,  Kinderman 
and  Ramage  [5]  use  Marsaglia's  [6]  algorithm,  as  modified  by  Ahrens  and 
Dieter  [1],  for  the  tail  of  the  normal  distribution.  In  the  context  of  the 
general  acceptance/rejection  algorithm,  Ahrens  and  Dieter  [2]  use  a normal 
majorizing  function  for  the  body  of  the  gamma  distribution  and  an  exponential 

majorizing  function  for  the  tails.  Schmeiser  and  Shalaby  [7]  use  a triangular 

majorizing  function  for  the  tails  of  the  beta  distribution. 

In  addition  to  use  as  part  of  other  algorithms,  tail  generators  are 

sometimes  used  directly.  The  most  straightforward  method  of  generating  from 
the  tail  of  a distribution  is  to  generate  a value  x from  the  complete  distri- 
bution and  use  x only  if  x lies  in  the  desired  tail,  but  this  method  is  inef- 
ficient for  small  tail  areas.  A second  method  is  available  when  the  inverse 
cumulative  distribution  function  transformation  is  used;  i.e.  x ■ F ^(u) 

where  u % U(0,1).  If  values  are  restricted  to  x >_  a,  the  transformation 
—2.  * * * 

x = F (u  + u(l  - u ))  is  used,  where  u = F(a).  Likewise  if  values  are 

-1  * 

restricted  to  x < a,  the  transformation  x = F (uu  ) is  used,  where  again 

* 

u = F(a).  Unfortunately,  such  common  distributions  as  the  normal,  gamma, 
and  beta  do  not  have  closed  form  inverse  cumulative  distribution  functions. 

A family  of  algorithms  is  developed  below  to  generate  variates  from 
distribution  tails.  The  algorithms  may  be  used  independently  or  as  part  of 
a complete  distribution  generator.  In  Section  2,  the  general  acceptance/ 
rejection  algorithm  is  specialized  to  use  an  exponential  majorizing  function 
to  generate  values  from  the  tails  of  distributions.  In  Section  3,  the 
normal,  gamma,  Weibull  and  beta  algorithms  are  given.  The  algorithms  are 
discussed  in  Section  4. 
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2.  THE  EXPONENTIAL  TAIL  ALGORITHM 


Let  f(x)  be  the  density  function  of  the  distribution  from  which  tail 
values  are  to  be  generated.  If  t(x)  is  a majorizing  function  of  f(x) 

(t(x)^  f(x)  for  all  x) , values  from  f(x)  may  be  generated  by  the  general 
acceptance/rejection  algorithm: 

1.  Generate  x having  density  r(x)  proportional  to  t(x). 

2.  Generate  v 'v  U (0,1). 

3.  If  v < f(x)/t(x),  deliver  x.  Otherwise  go  to  step  1. 

In  this  section,  this  algorithm  is  specialized  to  generate  only  tail  values. 

Given  an  appropriate  density  f(x),  the  figure  illustrates  the  algorithm 
developed  here  for  the  right  tail.  The  left  tail  algorithm,  discussed  later, 
is  quite  similar. 


Figure  About  Here 


To  ensure  the  algorithm  is  valid,  some  assumptions  on  f(x)  are 
necessary.  Assume  that  f(x)  is  convex  for  x > a,  which  is  true  of  "bell- 
shaped" distributions  if  the  point  a is  outside  the  points  of  inflection. 
Further  assume  the  tails  of  f(x)  are  lighter  than  the  exponential  tail. 
This  is  true  for  the  normal,  beta,  gamma,  and  Weibull  distributions,  but 
is  not  true  for  the  lognormal  or  Cauchy. 

An  exponential  majorizing  function  t(x)  ■ f(a)  exp[-A(x  - a)]  is  made 
tangent  to  f(x)  at  the  point  a by  setting  A » -f’(a)/f(a).  Since  t(x)  is 
proportional  to  the  exponential  density,  r(x)  ■ Ae  , x may  be 

generated  in  step  1 using  the  inverse  transformation 
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Commonly  step  3 is  implemented  using  logarithms,  so  here  step  3 
becomes 

If  In  v _<  J£n[f  (x)/t(x)  ] 

or 

Hn  v <_  2n(f (x)/f (a))  + X(x  - a). 

For  many  distributions  f(x)/f(a)  simplifies  nicely.  However  step  3 still 
requires  at  least  one,  and  commonly  more,  exponential  level  operations. 

Therefore  step  3 is  modified  to  include  a preliminary  comparison  of  v with 
b(x)/t(x)  before  comparing  v with  f(x)/t(x),  where  b(x)  3 f (a) [X(a-x)+l] 
is  the  linear  tangent  to  t(x)  and  f(x)  at  a.  Since  both  t(x)  and  f(x)  are 
convex  for  x _>  a,  b(x)  _<  t(x)  and  b(x)  f (x)  for  all  x a.  After  cancel- 
ling f (a) , the  preliminary  comparison  is 

v <_  (X[a  - x]  + 1) / exp [-X(x  -a)],  (1) 

which  simplifies  to 

v _<  [X(a  - x)  + l]/u 

after  substituting  x = a - £n(u)/Xinto  equation  (1). 

Allowing  the  steps  to  be  renumbered,  the  exponential  tail  algorithm  for 
a convex  density  f(x)  over  the  interval  x > a is  as  follows: 

1.  Generate  u,  v ^ U(0,1) 

2.  x * a - 2,n(u)/X 

3.  If  v < [X(a  - x) + l]/u,  deliver  x. 

4.  If  In  v <_  in[£ (x)/f (a)]  + X(x  - a),  deliver  x.  Otherwise  go  to  step  1. 
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The  algorithm  coding  is  unchanged  for  generating  values  from  the  left 
tail.  Assume  x <_  a and  f(x)  is  convex  over  this  region.  Then  X - -f'(a)/f(a) 
will  now  be  negative.  This  change  in  sign  for  X correctly  changes  each  step 
of  the  algorithm.  Thus  negative  values  of  X yield  the  left  tail  and  positive 
values  of  X yield  the  right  tail  automatically. 

A problem  does  arise  when  a negative  value  of  x is  generated  in  step  1 
for  distributions  where  f(x)  is  not  defined  for  negative  x.  In  this  case 
step  lb  can  be  added: 

lb.  u = e3^  + u(l  - e3^) 
aX 

where  e is  a constant  evaluated  only  one  time.  Equivalently,  step  2b  may  be 
added : 

2b.  If  x < 0,  go  to  step  1. 

| 

For  either  the  left  or  right  tail,  the  area  under  t(x)  is  of  interest 
when  generating  values  from  the  entire  distribution.  The  probability  of 

’ 

needing  a tail  value  for  a particular  variate  is  the  ratio  of  the  area  under 
the  tail  majorizing  function  to  the  area  under  the  majorizing  function  for 
the  entire  distribution.  For  the  exponential  tail  algorithm  discussed  here, 
tail  area  is  |f(a)/X|  for  either  tail. 
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3.  SPECIALIZATION  TO  SOME  COMMON  DISTRIBUTION '» 

Specialization  to  the  normal,  gamma,  Weibull,  and  beta  distributions 
is  now  straightforward.  Let  K(x)  denote  Jln[f  (x)/f  (a)  ] . Scaling  constants 
in  f(x)  are  of  no  consequence  here  and  are  ignored. 

2 

Normal  f(x)  ■ exp(-x  /2)  - °°  < x < °° 

I a | > 1 
X * a 

K(x)  = -(x  - a)2/a. 

Gamma  f(x)  =*  x°  ^"exp(-x)  0 x,  a 1 

If  a 2,  there  is  no  left  tail.  Otherwise 

0 < a < (a  - 1)  - (a  - 1)1/2  or  (a  - 1)  + (a  - 1)1/2  < a 
A = 1 - [(a  - 1) /a] 

K(x)  * (a  - l)£n(x/a)  + a - x 

Weibull  f(x)  = xC  ^exp[-xC]  0 <_  x,  c 1 

If  c £ 2 there  is  no  left  tail.  Otherwise 
0 < a < { [ (c  - 1)  - (c  - l)1/2]/c}1/c  or 

{[(c  - 1)  + (c  - l)1/2]/c}1/c  < a 

A » cac  1 - (c  - l)/a 

K(x)  ■ (c  - l)[Zn(x)  - &n(a)]  + aC  - exp[c  Zn(x)] 

Beta  f(x)  * xP(l  - x)^  0 <_  x _<  1,  P>0,  Q>0 

If  P + Q < 1,  there  are  no  tails. 

If  P + Q > 1,  then  let  D - (PQ/(P  + Q - 1))1/2 

If  P - D < 0,  there  is  no  left  tail,  and  if  Q - D <_  0, 

there  is  no  right  tail. 

Otherwise,  0 < a <_  (P  - D)/(P  + Q)  or  (P  + D)/(P  + Q)  <_  a < 1. 
A - Q/(l  - a)  - P/a 

K(x)  = P 4n(x/a)  + Q l n((l  - x)/(l  - a)) 
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4.  ANALYSIS  OF  THE  EXPONENTIAL  TAIL  ALGORITHM 


A result  which  simplifies  the  analysis  of  the  algorithms  is  that  the 
linear  comparison  of  step  3 will  be  true  with  probability  .5  for  any  f(x) 
and  for  any  a.  This  is  easy  to  verify,  since  the  area  under  b(x)is  |f(a)/(2A)| 
and  the  area  under  t(x)  is  ] f ( a) /X | . 

Since  the  exponential  level  operations  dominate  the  computation  times 
on  most  computers,  crude  algorithm  comparisons  can  be  made  by  counting  the 
expected  number  of  such  operations.  Thus  the  exponential  tail  algorithm 
requires  n * 1 + 1 + (m/2)  expected  logarithms  where  m is  the  number  of 
logarithms  required  to  evaluate  K(x) . The  four  distributions  discussed  in 
Section  3 result  in  m » 0,  1,  2,  and  2 respectively.  (The  two  exponential 
operations  in  the  Weibull  algorithm  are  £n(x)  and  exp(-)).  The  expected  total 
operations  are  then  n * 2,  2.5,  3,  and  3,  respectively. 

The  other  factor  in  determining  execution  time  is  the  expected  number 

00 

of  iterations  per  variate  generated,  which  is  [f(a)/A ] //  f(x)dx  for  the 
right  tail  and  [-f(a)/A ]//^  f(x)dx  for  the  left  tail.  For  the  four  distri- 
butions considered  in  Section  3 the  expected  number  of  iterations  is 
relatively  constant  for  all  values  of  a,  due  to  all  four  distributions  having 
tails  not  unlike  the  exponential. 

It  is  interesting  to  compare  Marsaglia's  [6]  normal  tail  generator  as 
modified  by  Ahrens  and  Dieter  [1]  to  the  algorithm  of  this  paper.  Comparison 
was  performed  on  the  SMU  CDC  Cyber  72  Computer  under  the  KRONOS  operating 
system  and  using  the  FTN  Fortran  compiler  with  its  internal  source  of  17(0,1) 
values  RANF.  Both  are  easy  to  implement,  with  Marsaglia's  algorithm  requiring 
seven  lines  of  code  and  the  algorithm  of  this  paper  requiring  six  lines. 
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The  algorithms  ran  with  almost  identical  times,  ranging  from  .45  milli- 
seconds per  variate  for  a - 1 to  .33  milliseconds  for  a » 5.5.  The  ratio 
of  uniform  dev  iates  method,  given  by  Kinderman  and  Monahan  [4],  yields  a 
Normal  tail  algorithm,  described  in  Kinderman  and  Monahan  [3],  which 
generates  variates  faster  but  which  requires  a one-time  initialization 
of  constants. 
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